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SUMMARY 


Five viscous transonic airfoil cases have been computed by two significantly different compu- 
tational fluid dynamics codes: an explicit finite-volume algorithm with multignd, and an implicit 
finite-difference approximate-factorization method with eigenvector 
are described in detail, and their performance on the test cases is compar . 
same grids, turbulence model, and computer to provide the truest test of the algo t . P 

proaches produce very similar results, which, for attached flows, also agree well with expenmenta 

results; however, the explicit code is considerably faster. 


INTRODUCTION 


The development of algorithms for the Navier-Stokes equations has proceeded along two distinct 
bu, intenwtnedpaths: exploit methods, which are iess compuut.ktnally intensive but ate reacted by 
stability considerations to small time-steps, and implicit methods, which trade a larger operation count 
to meLer smbUto. Jameson e. al. pioneered the use of an explicit modified Runge-Ku.ta method for 
the Euler equations (ref. 1). With the enhancements of multigrid, implicit residual smoot tng, an 
local time-stepping, the method was practical for application to the Navier-Stokes equations, na o- 
gously Beam and Warming developed an implicit approximate-factorization scheme which was usefu 
fwreducing die work of a two- 01 three-dimensional implicit operator to a senes of one-d,mens.onal 
operators (rtf. 2). Enhanced by local time-stepping, eigenvector diagonal, zation, ^ ^2 
which applies some of the ideas of multigrid, this algorithm also proved applicable to the Navier-Stokes 

equations. 

Both methods have been rigorously tested against experimental data for a variety of flow condi- 
tions and both are in widespread use today in research and practical applications. The purpose of this 
study is to compare the methods with regard to their capability to compute transonic viscous a 
flows Computations were performed using two codes: FLOMG, which is based on Jameson s mulu- 
grid algorithm, and ARC2D, which uses the Beam-Warming implicit method. Since compumtiomd 
results can be sensitive to various features of the grids, both codes used the same^gnds^hree of which 
were generated for each airfoil. To provide the most controlled test of the methods, both codes als 
used the same turbulence model and computer. Five cases, with conditions ranging from subcntical 
attached flow to strong shock-induced separation, were computed. Results include plots of Pressure 
coefficient skin-friction coefficient, boundary-layer profiles, and displacement thickness, as we 

Efficients and convergence times. Comparisons with experimental data are made 

when data are available. 
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NAVIER-STOKES EQUATIONS 


The nondimensional, strong-conservation-law form of the two-dimensional Navier-Stokes equa- 
tions in Cartesian coordinates is 


d t Q + d x E + d v F = Re~\d x E v + d v F v ) 
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and Re is the Reynolds number, p is the density, u and v are the Cartesian velocities, and e is the total 
energy, and where 

fix = P ( 4 Uj 2t>y)/3 

T xy = p(u„ + u z ) 

t vv = p(- 2 u 2 + 4 v v )/3 ( 2 c ) 

/t = t it xx + + pPr -, ( 7 — l) _ 1 dz 0 2 

Sto = + vr tfV + pPr -1 (7 - l) - 1 d„a 2 

Here Fr is the Prandtl number, and a is the speed of sound (a 2 = 7 p/p for ideal fluids). 

Pressure is related to the conservative flow variables Q by the equation of state. 


P = (7~ 1) 


e - jp(u 2 + v 2 ) 


(3) 


where 7 is the ratio of specific heats, generally taken as 1.4, and p is the dynamic viscosity, which is 
made up of the molecular viscosity plus a computed turbulent eddy viscosity. 

The choice of nondimensional parameters is arbitrary. Here we have chosen to scale the variables 
as 

~P~ti_v_ e 

p= — , u= — , v = — , e = T ( 4 a) 

Poo Ooo ®oo Poo®oo 

where 00 refers to free-stream quantities. Assuming a reference length l (usually taken as some char- 
acteristic physical length such as chord of an airfoil), time t scales as t - ta/l. The viscous coefficients 
scale as 



be dropped for simplicity. 

Coordinate Transformation 

For a finite-difference solution on a body-fitted grid, the equations are transformed “ described 
in references 3 and 4. The Navier-Stokes equations written in generalized curvtl.near coordtnat (£, 

° ale d ,Q + d ( E + 3 ,F = Re'' l d ( E> + 3„FJ < 5 > 
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with J-' = <*{¥, -*,»(> representing the metric Jacobian and U - 

v representing the contravariant velocities. The viscous flux terms are E, = J ( i.B, + V,) 
F„ = J~ l (t) x E v + 1 TyF v ). 

The stress terms, such as r„, are also transformed in terms of the £- and ^-derivatives where 

T xx - + Vx u v^ — 2(^ + 

T xy = n(£yUt + Vy u V + + *1 xV ^ 

Tyy = /j[-2(£xU£ + V x u v ) + 4 + T) v Vt,)]/3 

/4 = uTu + + pPr~ x ( 7 — 1) l ($x^£a + rj x d v a ) 

g 4 = u T xy + i>t„ v + /iPr -1 (Tf - l) -1 (£yd£° 2 + %3r,a 2 ) 


(7) 


Thin-Layer Approximation 

In high-Reynolds-number viscous flows, the effects of viscosity are concentrated near rigid bound- 
aries a„d1n wake regions. Typical grids have finer grid spacing in directions nearly notmaltothe 
surfaces and coarser grid spacing along the surface. On such grids, the Viacom t terms 
derivatives along the body will not be resolved, and in most cases for attached and ™ d y p , 
flows these terms are negligible. The terms associated with the normal direction will be resolved for 

sufficiently fine-grid spacing, and these are substantial terms. 

The thin-layer approximation is similar in philosophy to, but not the same as, the boundary-layer 
theory for which appropriate scaling arguments show that streamwise components of the viscous terms 
neglected relative to the normal terms. In the thin-layer approximation, the normal momentum 
equation is solved and pressure can vary through the boundary layer. 
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Applying die thin-layer approximation to equations (5)-(7), in which all the viscous terms associ- 
ated with £ -derivatives are neglected, the following system of equations is obtained: 

d T Q + d ( E + d v F = Re-' d v S ( 8 ) 


where 
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TUrbulence Model 

The algebraic mixing length model of Baldwin and Lomax (ref. 5) is included to approximate 
the effects of turbulence. It is a two-layer model in which the inner layer is governed by the Prandtl 
mixing length with Van Driest damping; the outer layer follows the Clauser approximation. Computed 
vorticity is used in defining the reference mixing length required for the outer layer. In order to improve 
the numerical compatibility of the Baldwin-Lomax model, two simple modi fi cations are made. First, in 
the Van Driest damping factor, the shear stress at the wall is replaced with the maximum local laminar 
shear stress. This generally prevents the numerical difficulties caused by vanishing shear stress at 
separation. Second, the constant C wk , which is used in the outer-layer viscosity formulation, is changed 
from 0.25 to 1.0. For some transonic flows, convergence may not be possible if this value is too low, 
a result of motion of the shock wave. 


NUMERICAL ALGORITHM FOR FLOMG RESULTS 


The FLOMG computations were performed using a Navier-Stokes code developed by Swanson 
and Turkel, which is based on the explicit multistage time-stepping schemes of references 1 and 6 . This 
class of schemes is currently in widespread use for solving the Euler equations of gas dynamics. In 
references 7-10, these schemes were extended to allow the solution of the compressible Navier-Stokes 
equations. Significant improvements in numerical efficiency were introduced in references 11-13. The 
basic ideas developed in two dimensions were extended to three dimensions in references 14-16. In 
the code of Swanson and Turkel, a cell-centered, finite-volume method is employed to obtain differ- 
ence approximations for the flow equations. Such a method provides flexibility in treating arbitrary 
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geometries and different grid topologies, since no special treatment is ' ‘^^""^‘"(Se 

i^^MorTdefi^iti^^r'suffi^iemly^^tw^i-^rtificLl dissipation terms with adaptive ^retring 

™„^™cSeml ."es These merits include local time-stepping (a precond.uon.ng 

r,:rr:e^^ 

viscous flow problems; they also work quite well on typ.cal meshes for mv.sctd flows. 


Spatial Discretization 


A finite-volume approach is applied to discretize 
domain is divided into quadrilateral cells, fixed in time, 
written in integral form as follows, 


the equations of motion. The computational 
and for each cell the governing equations are 


— f [ Qdxdy + [ ( Edy - Fdx) = Re 1 (E v dy - F v dx) 
dtJJn Jsa Jdn 


where Q is a generic cell and dft its boundary. 

A cell-centered type of discretization is used, and the line integral of equation ( 1 0) is approximated 
with the midpoint rulel^Thus, taking Q,y as a cell-averaged solution vector, equation (10) can be written 

in semidiscrete form as J (11) 


dt 




where £ is a discrete operator defined by 


C = Cc + Fn ~~ Fad 


with the subscripts C D, and AD referring to convection, diffusion, and artificial dissipation, re- 
r^hvely Ae conve'c.ive fluxes a. the cell faces are obtained by an averaging process. Moreover, 

summing over the cell faces, 

CcQ<j-'E*- § ‘ <12) 


1 - 1 


with the flux tensor, which is associated with convection, given by 

Fi = Eie x + Fie y , 


and for each cell face l, the directed area S, is expressed as S, = ( A y), & — ( A *) , e,, where the proper 
signs of (A x)i and ( Ay) , produce an outward normal to the cell face. The augmented convecuve flux 

tensor is evaluated as 




-V- + Q + V + ) t + Pi 


( 13 ) 
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Figure 1 . Finite- volume discretization. 

~ | ® (Pavg)i^x (Pavg) C(P^)ot;j)j j 

(Pa„,)i= j(p- + P*) t , ((PVW,= J(P~V ~ + P + V + )i 

and the minus and plus superscripts mean that one quantity is determined at the center of the cell of 
interest and the other at the center of the adjacent cell on edge l, and V is the velocity vector. The 
spatial derivatives necessary to compute the viscous terms are evaluated by means of Green’s theorem. 
For example, consider the arbitrary cell ( ABCD ) in figure 1. The contributions u x and u v to the 

viscous flux across each cell face (i.e., BC) are approximated with their mean values using 



where Q' is an auxiliary cell {A! B’C’D' for face BC in fig. 1). The values of u at B and C required 
for the cell-boundary integrals in equation (14) are obtained by a simple average of the values of u at 
the four surrounding cells. In a similar manner, v x and v v are computed. Additional details for finite- 
volume treatment of viscous stresses and heat-conduction terms are given in reference 7, which also 
gives a thin-layer formulation, and in reference 9. 
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Artificial Dissipation 


For inviscid flows the basic finite-volume scheme just described contains no dissipative terms. 
In order to prevent odd-even point decoupling and oscillations near shock waves or stagnation points, 
artificial dissipation terms are added to the governing discrete equations. Moreover, the introduction 
of appropriate dissipation in the vicinity of shocks allows an entropy condition to be satisfied, and thus 
guarantees the uniqueness of weak solutions. For viscous flows, dissipative properties are present ow- 
ing to diffusive terms; however, due to nonlinear effects the physical dissipation may not be sufficient 
to guarantee stability, especially in the case of the highly stretched meshes generally used to resolve 
the steep gradients in shear layers. Thus, to maintain the stability and robustness of the numerical 
procedure, artificial dissipation is also included in viscous regions. 

The artificial dissipation model used in the scheme is basically the one developed by Jameson et al. 
(ref. 1). This nonlinear model is a blending of second and fourth differences. The quantity LadQ\J in 
equation ( 1 1 ) is defined as 


C'AoQij - ( — D* + D* — D v ) Qij 

(15) 

where (£, 17 ) are arbitrary curvilinear coordinates. 


D\Q,j = v £ < Vw ’tfijWiJ 

(16) 

D(Qij = E !*V AfVfAfQiJ 

(17) 


and where i,j are indices (for a cell center) associated with the £- and rj-directions, and are 

forward and backward difference operators in the £ -direction. Following references 11 and 12, the 
variable scaling factor is defined as 


where 

(^£ )ij = <M r > 

(19) 

= 1 

and where r is the ratio A and 316 the scaled spectral radii of the flux Jacobian matrices 

(associated with the £- and rj-directions) for the Euler equations, and the exponent < is generally taken 
to be 2/3. The spectral radii for the £- and 77 -directions are given by 


= \uy v - vx v \ + ayfy* + x\ 


= |vx£ - uj/{| + ayjx\ + yj 
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and a is the speed of sound. The coefficients e (2) and e (4) use the pressure as a sensor for shocks and 
stagnation points, and they are defined as 


e<?, , = max(ui-\j,vij,vn.\j,ui+2j) 

y J 


( 22 ) 
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ij ~ 


Pi-Vj ~ 2 Py + Pi+lj 
Pi-IJ + 2 PiJ + Pi+ 1 J 



max[0,(* (4) - £ f (2 \ j)] 


(23) 

(24) 


where typical values for the constants « (2 ^ and are in the ranges 1/4 to 1/2 and 1/64 to 1/32, 
respectively. For the normal direction ( rj) , the dissipation contributions are defined in a similar way, 
except 

(25) 


~ ) ( ^*))i 


V'*J 


The treatment of the artificial dissipation must be modified at the boundaries of the physical do- 
main. In the case of the fourth-difference dissipation, the standard five-point difference stencil must be 
replaced at the first two interior mesh cells. This means that one-sided or one-sided biased stencils are 
used at these cells. The dissipative character of the artificial terms is important because it influences 
both stability and accuracy. For example, if the dissipation is too large at a solid boundary, an artificial 
boundary layer is created in an inviscid flow, and the effective Reynolds number for a viscous flow is 
altered. To improve accuracy at the wall boundaries of viscous flows, where gradients are steep be- 
cause of physical boundary layers, the usual fourth-difference stencils are changed in this dissipation 
model. 

Let the total dissipation for a mesh cell, in the direction represented by the index j, be denoted by 
dj. For simplicity, assume that Xe <4) = 1 . Then, 

(26) 


where the dissipative flux 


df J+ . = (AQ) /+ | - 2(A Q), + . + (A Q),._. (27) 

and thus 

dj = (AQ) ; ,3 - 3(AQ) ;> , + 3(AQ) ; ._. - (AQ) ; _, (28) 

with the index i for A Q suppressed for convenience. Consider the first two interior cells adjacent to a 
solid boundary, as depicted in figure 2. If 


(AQ). =(AQ)b =(AQ)| 


then equation (28) gives 


dz = Qa — 2Q3 + Q2 
di =Q 5 -4Q 4 +5Q3-2Q 2 


(29) 
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(30) 

(31) 




Figure 2. Boundary point dissipation. 


These boundary stencils are fairly standard ones, and they result in a nonpositive definite dissipation 
matrix for the system of difference equations (ref. 17). An alternative form, which has reduced the 
sensitivity to solid surface normal mesh spacing for turbulent flow calculations without compromising 
stability or convergence, is given by 

(AQ)i = 2(AQ)3 - (AQ)s (32) 


and 

d .2 = Q4 — 3 Q 3 + 3Q2 — Q 1 
di = Q 5 — 4 Q 4 + 6 Q 3 — 4 Q 2 + Qi 


(33) 

(34) 


Time-Stepping Scheme 


The system of ordinary differential equations of equation (11) is advanced in time toward the 
steady-state solution with a five-stage Runge— Kutta scheme. This scheme is second-order accurate in 
time. At the (m + l)sf stage, 




= QS? -a m+ ,^£[ C c Q\ 
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where = Q", is the discrete solution at time-level n; is the solution at the new 

time-levef, n+ 1; a m+ i are the coefficients of the scheme; At is the time-step; £2 is the area of the 
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mesh cell; and 7 mn are the weighting factors of the artificial dissipation. The coefficients <* ra+ i are 
determined such that the scheme has the largest possible hyperbolic stability limit. An appropriate set 
of coefficients is given by 


113 1 

= 4’ “ 2= 6’ “ 3= 8’ “ 4 = r 


(*5 = 1 


(36) 


This scheme also exhibits good high-frequency damping behavior, which is crucial for a rapidly conver- 
gent multigrid method (to be described in a subsequent section). In order to establish a good parabolic 
stability, the artificial dissipation terms are evaluated on the first, third, and fifth stages. The weighting 
factors 7 TOn must satisfy the condition 

07) 

They are defined as follows: 


Too = 1, 

7io = 1 , Til = 

720 = ( 1 - 73 ) . 721 = o , 722 = 73 , (38) 

730 = ( 1 - 73 ) » 731 = 0 , 732 = 73 . 733 = 0 , 

740 = (1 -73)0 -Ts>> 741= 0, 742 = 73(1 7s ) > 743 = 0, 744 = 75 

where 73 = 0 .56 , and 75 = 0 .44 . As indicated in equation (35), the physical viscous terms are 
computed only on the first stage and are frozen for the remaining stages. The single evaluation appears 
to have no significant effect on the stability of the scheme and allows a reduction in computational 
effort. 


Convergence Acceleration Techniques 

Three methods are employed to accelerate convergence of the basic explicit time-stepping scheme. 
These techniques are: (1) local time-stepping; (2) residual smoothing; and (3) multigrid. With local 
time-stepping, the solution at each mesh point is advanced at the maximum At allowed by stabil- 
ity. Both convection and diffusion limits are included in the At computation (see ref. 18). Implicit 
smoothing of the residuals is used to extend the stability range of the basic time-stepping scheme. For 
two-dimensional problems, the residual smoothing can be applied in the form 

(1 -ftV^Xl (39) 

where the residual is defined by 

■Rij' + CdQ^ - AD fm) ) , m = 1,5 (40) 

and is computed in the Runge-Kutta stage m, AD (m) is the total artificial dissipation at stage m, and 
7£ (m) is the final residual at stage m after the sequence of smoothings in the (- and rj-directions. The 
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coefficients /% and (3 V are variable and are functions of the spectral radii X £ and X n . Based on the 
unpublished work of Swanson, they can be written as follows: 


/% = max 



JL 1 —V _ i 

N* 1 + 



= max 



- — !-tY - 1 

N- l + *r 3) 



(41) 


where the ratio r* = X,/X £ , and the quantity N/N* is the ratioof the Courant-Friedrichs-Lewy number 
of the smoothed scheme to that of the basic explicit scheme (usually having a value of 2). From a 
linear stability analysis, the scheme with these coefficients is stable for all mesh-cell aspect ratios when 
the parameter V> « 0 125 and N/N* is sufficiently large. The practical limitation on the Courant 
number is due to the requirement for effective high-frequency damping. For large N/N*, the high- 
frequency damping of the scheme vanishes. An alternative form for these coefficients was introduced 
by Martinelli and Jameson (ref. 1 2). Similar performance of the two forms has been observed for typical 
viscous flow meshes. However, significant improvement in convergence rate has been obtained with 
standard inviscid meshes using the form of equation (41), 

The multigrid method used with the Runge-Kutta time-stepping algorithm is based on the work 
of Jameson (ref. 6). The basic idea of the method is to use coarser grids to speed up the propagation 
of the fine-grid corrections. Coarser meshes are obtained by eliminating every other mesh line in each 
coordinate direction. On the auxiliary meshes, the solution is initialized as 


n (0) _ 

Qlh " n 2h 


(42) 


where the subscript refers to the mesh-spacing value, the sum is over the four fine-grid cells that com- 
pose the 2 h grid cell, and, again, Q is a cell volume. This rule conserves mass, momentum and 
energy. On a coarse grid, a forcing function P is added to the governing discrete equations in order to 
impose the fine-grid approximation. After the initialization of the coarse-grid solution, this function is 


computed as follows: 


(43) 


where Rh(Qh) = ChQh- Then, the time-stepping scheme on the (m + l)st stage becomes 


n (m+\) 




(44) 


We can also define a new value f?* for the residual as 

R\ h = RihiQzh) + Pih ( 45) 

This value can be collected, the solution Q 2/1 restricted to the next coarser grid, and the process re- 
peated. The corrections computed on a coarse grid are transferred back to a finer grid with bilinear 
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interpolation. A fixed W - type cycle is used to execute the multigrid strategy (see ref. 18 for a diagram 
of a cycle). In order to make this strategy effective for a wide range of flow conditions, the resultant 
coarse-grid corrections are smoothed before they are passed to the finest mesh. The factored scheme 
of equation (39) with constant coefficients ($ = ^ = 0 .2 ) is used for this smoothing. Also, the ap- 
plication of a full multi grid ( FMG ) method provides a well-conditioned starting solution for the finest 
mesh being considered. Additional details of the multigrid procedure are discussed in reference 18. 


Discrete Boundary Conditions 


At a solid boundary, a row of auxiliary cells is created exterior to the domain of the flow. The 
no-slip condition is imposed by treating the Cartesian velocity components as antisymmetric functions 
with respect to the solid surface. Thus, 

U,,l = — Uj t 2 


(46) 


V,.l = —Vi, 2 


where the indices (i, 1) and (i, 2) refer to the centers of the auxiliary and the first interior cells, re- 
spectively. The surface values of pressure p and temperature T are computed using the reduced normal 
momentum and energy equations. 




(47) 


where rj is the coordinate normal to the surface. 


At the outer boundary of “C” meshes, the boundary points are treated in essentially the same way 
as described in the discussion of boundary conditions in the next section. Simple extrapolation is used 
to compute the flow quantities at the downstream boundary of the C-type mesh topology. 


NUMERICAL ALGORITHM FOR ARC2D RESULTS 


The ARC2D computations were performed using a code based on the implicit approximate factor- 
ization algorithm of Beam and Warming (ref. 2). The code was originally developed by Steger (ref. 19) 
and has been steadily improved and modified (refs. 20-22). The algorithm is an implicit approximate- 
factorization finite-difference scheme which is first-order accurate in time. Local time-linearizations 
are applied to the nonlinear terms and an approximate factorization of the two-dimensional implicit 
operator is used to produce one-dimensional operators. The spatial derivative terms are approximated 
with second-order central differences. A diagonal form of the algorithm that produces a computa- 
tionally efficient modification of the standard algorithm in which the diagonalization results in scalar 
pentadiagonal operators in place of the block operators is also used. Explicit and implicit artificial 
dissipation terms are added to achieve nonlinear stability. Spatially variable time-steps and mesh se- 
quencing are used to accelerate convergence. 
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Implicit Time-Differencing and Local Time-Linearization 


Applying the first-order Euler implicit time-differencing scheme to equation (8) results in 

Q*' -Q*+ h(E?' + Ff' - Re-'S?') = 0 (48) 

with h = At. The nonlinear terms are linearized in time about Q n by a Taylor series such that 


E** 1 = E n + A n AQ n + 0(h 2 ) 
_ pn + § n AQ n + 0{h 2 ) 


=Re~ x [$* + J~ x M n AQ n ) 


+ 0(h 2 ) 


(49) 


where A = dE/dQ, B - dF/dQ, and M = dS/dQ are the flux Jacobians and AQ n is 0(h). The 
linearizations are second-order accurate. 


The Jacobian matrices are A or B = 

K t « I ^ 

—ud + K x cf > 2 K t + 0 - (7 — 2 )k x u K y u — (7 — 1)k z v (7— 1)«i (5Q) 

—v 6 + k v 4> 2 k x v — (7 — l)« v u + 6 — (7 — 2 )k v v (7— 1)k v 

. 6 [(f > 2 — ui ] « a O! - (7- 1)U0 -(7- \)v 6 7 - 

with 01 = 7(e/p) - (f> 2 , 6 = k x u + k v v, (j > 2 = 3/7 - l)(u 2 + v 2 ), and k = ( or 77 for A or B, 
respectively. 


The viscous flux Jacobian is 

- 0 

M= m21 
m 3 i 

_m 4 i 


0 0 
a\d v (p~ x ) a 2 d^(p~ x ) 
<*2 d„(p~ x ) «3 d v (p~ ] ) 


77142 


17143 


0 
0 
0 

77144 J 


(51a) 


77X21 = - a\d v (u/p) - a 2 dr,(v/p) 

77X31 = — Ol 2 d v (u/p) — Ct3d v (v/p) 

77X41 =014 dr, [~ (e/p 2 ) + ( U 2 + V 2 )/p] 

- a\d n (u 2 /p) - 2 CK 2 d v ( uv/ p) 

- OiT,dr,(v 2 Ip) 

77X42 = - a 4 9„(u/p) - 77X21 
77X43 = - 0i4d v (v/p) - 77131 
77X44 =Oi4dr,(p~ X ) 

0t\ =/z[ ( 4 /3) rjx 2 + i) y 2 ], oti = (p/^VxVv 

a 3 =p[r)x 2 + (4/3) 77 y 2 ] , a 4 = lpPr~\r} x 2 + T} y 2 ) 


(516) 
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Applying equations (49) to equation (48) and combining the AQ n terms produces the unfactored 
block algorithm in delta form. 


[/ + hd^A n + hd„B n - Re~ l hJ-'d v M^Q n = 
-h (dtE* + 3„F n - Re~ x d v S n ^) 


(52) 


Space Differencing 

The next step is to take the continuous differential operators dt and d v and approximate them 
with finite-difference operators on a discrete mesh. Introducing a grid of mesh points (;, k) , variables 
are defined at mesh points as u ;> * = u(;'A£, fcA rj) . The grid spacing in the computational domain is 
chosen to be unity so that A£, A tj = 1 . Second-order central-difference operators are used, where, for 
example, 

= (u ;+ i ,fc - /2 and 8 v Uj, k = (uy.it+i - / 2 < 53 °) 

For the viscous derivatives, the terms take the form 

d n (a jik d v /3j, k ) (53 b) 

which is differenced in the compact three-point form as 

[(a,>i + «>,*) (PiMi - Pit) - (“>,* + (#.* - Pit- 1)] / 2 (53c) 


Approximate Factorization 


The implicit (left-hand) side of equation (52) can be written as 

(i + h8 t A n + h8 v B n - hRe~ x 8r,J~ x M n ) A Q n = 
(/ + h8 ( A”) (/ + h6 n B n - hRe~ x 8 V J ~ x M n ) AQ" 
-h 2 8 ( A n 8 v B n AQ n + h 2 Re-'8 ( A n 8 v J-'M n AQ a 


The cross term is second-order accurate, since A Q n is 0(h). It can therefore be neglected without 
degrading the time-accuracy of any first- or second-order scheme. 

The resulting factored form of the algorithm is 


(/+ (/ + h8 v B n - hRe-'S^J-' M B ) AQ n = 

-h + 6 r) F n - Re~ x 6„S n ) 


(55) 
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with two implicit operators, each of which is block tridiagonal. The solution algorithm now consists of 
two one-dimensional sweeps, one in the £- and one in the rj-direction. Each sweep requires the solution 
of a linear block- tridiagonal system by block LU (lower-upper) decomposition. 

Diagonal Form of Implicit Algorithm 

One way to improve the efficiency of this method is to reduce the operation count by introducing a 
diagonalization of the blocks in the implicit operators, as developed by Pulliam and Chaussee (ref. 23). 
The eigensystems of the flux Jacobians A and B are used in this construction. Initially, the discussion 
is restricted to the Euler equations; the application to Navier-Stokes is considered later. 

Each of the flux Jacobians A and B has real eigenvalues and a complete set of eigenvectors. 
Therefore, the Jacobian matrices can be diagonalized (see ref. 24) 

A i = T{'AT i and A, = T~ ] BT V 

with Ti the matrix whose columns are the eigenvectors of A and T v the corresponding eigenvector 
matrix for B. 

Replacing A and B in equation (52) with their eigensystem decompositions produces 

[r^ + fcfi^Aerf 1 )] [t„t~ 1 + hs v (T 1l A.r- 1 )] a q* 

= the explicit right-hand side of equation (52) = R n ( 56) 

At this point, equation (56) and the inviscid form of equation (55) are exactly equivalent. A 
modified form of equation (56) can be obtained by factoring the and T n eigenvector matrices outside 
the spatial derivative terms 8$ and 8 n . The eigenvector matrices are functions of £ and rj and, therefore, 
this modification reduces the time-accuracy to at most first-order. The resulting equations are 

Tt(l+h8 ( A*) N (I + h 8 V \) T~ i AQ n = R n (57) 


where N -T^ x T n . 

The diagonal algorithm as presented in equation (57) is only rigorously valid for the Euler equa- 
tions. The implicit linearization of the viscous flux S n in the implicit operator for the redirection was 
neglected, because the viscous flux Jacobian M n is not simultaneously diagonalizable with the inviscid 
flux Jacobian B n . For stability and robustness, it is desirable to include the viscous term on the implicit 
side, but retaining the block operator in the redirection requires more computation time. Adding an 
approximate viscous eigenvalue. 


A„(r j) = p/iHe 'j '(fJi + f? y) 


(58) 


to the implicit operator has proved to be an appropriate compromise between accuracy and speed. 
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The explicit side of the diagonal algorithm is exactly the same as in the original algorithm, equa- 
tion (52). The modifications are restricted to the implicit side and so, if the diagonal algorithm con- 
verges, the steady- state solution will be identical to one obtained with the unmodified algorithm. The 
diagonal algorithm reduces the block-tridiagonal inversion to 4 x 4 matrix multiplies and scalar tridi- 
agonal inversions. The overall savings in computational work can be as high as 40%. 

Nonlinear Artificial Dissipation Model 


A mixed second- and fourth-difference dissipation model with appropriate coefficients produces 
a scheme with good shock-capturing capabilities. Jameson et al. have employed a dissipative model 
of such a form in which second- and fourth-difference dissipation are combined (ref. 1). The model, 
written in our notation, is 


(l/V + 1 ,* JjA t k + Jj,k ) C j,k 


(59) 


with 


=/t 2 Af max(/i ;+ i i *,/x ; _i,fc) 


_ jgriM ~ + PjzlA 1 ( 60) 

^' k |P;'+l,Jk + ip/.Jfe + Pj-l,k\ 

= max ( 0 , «4 A t - ) 

where for these calculations values of the constants are k 2 = 1 and K 4 = 1/100. Similar terms are 
used in the Tj-direction. The term ^Pj t k Is a spectral radius scaling factor and is defined as 


t/»i =|(/| + +e v (61) 
V'y =1^1 + a -]/nl + vj 


which is the sum of the spectral radii of A and B. 

The first term of equation (58) is a second-difference dissipation term with an extra pressure gradi- 
ent coefficient to increase its value near shocks. The second term is a fourth-difference term where the 
logic to compute switches it off when the second-difference nonlinear coefficient is larger then the 
constant fourth-difference coefficient. This occurs near a shock. This model is added to the right-hand 
explicit side of the algorithm. 

The implicit dissipation used with equation (58) is the linearization of the model, treating the 
pressure coefficient n and the spectral radius i/> as space-varying functions but ignoring their func- 
tional dependency on Q. Then the dissipation is linear in Qj t k *tnd is added to the diagonal algorithm, 
necessitating scalar pentadiagonal solvers. This produces the most efficient, stable, and convergent 
form of the implicit algorithm. 

Near computational boundaries, the fourth-difference dissipation is modified so as to eliminate 
the five-point stencil. A one-side biased second-difference term is used instead. 
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Convergence Acceleration 


In solving these steady-state problems, two techniques were used to accelerate the algorithm, the 
goal being to eliminate the transient as quickly as possible. Note that for the delta form of the algorithm 
(either factored or unfactored) the steady-state solution is independent of the time-step, h. Therefore, 
the time-step path to the steady-state does not affect the final solution and we can use time-step se- 
quences or spatially variable time-steps to accelerate convergence. The particular form of spatially 
variable time-step used here is to replace h with 



with A f 0 chosen to be 0{ 1) . The use of the local time-step can improve convergence time by a factor 
of 2 to 3. 

Another way to accelerate convergence to a steady state is to obtain a good initial guess for a 
fine mesh by first iterating on a sequence of coarse grids and then interpolating the solution up to the 
next refined grid. A coarsened grid is cut from each previous grid by halving the number of points 
in the ^-direction and by regenerating a new redistribution using fewer points. The redistribution is 
regenerated because in viscous flows the spacing near the wall would be too coarse if the halving 
procedure were used. A small number of iterations (a few hundred) are carried out on each coarsened 
grid, then the approximate solution is interpolated onto a more refined grid. The finest grid is then 
iterated to convergence. This can reduce the practical convergence time by as much as a factor of 2. 


Boundary Conditions 

At a rigid-body surface, the no-slip condition is satisfied with the Cartesian velocities u , v set to 
zero. The pressure is obtained through extrapolation such that dP / dr] = 0. A zeroth-order extrap- 
olation is used to compute density. At subsonic free-stream conditions, the outer boundary condition 
is based on locally one-dimensional Riemann invariants. At a boundary, local normal and tangential 
velocity components are computed. The normal component is the scaled contravariant velocity, 

y fl = 4SL (63) 

and the tangential component is 

Vt _ M ~ (64) 

^/(Vi + T Jy) 

The locally one-dimensional Riemann invariants are 

Rj = V n - 2 o /('7 - 1) and R 2 = V„ + 2 0/(7 - 1) (65) 

Two other equations are needed so that four unknowns (the four flow variables) can be calculated; Vt 
and a quantity S = p/p 1 , which has the same functional dependence as entropy, are used. For an inflow 
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boundary, V n < 0 and then three quantities can be specified. The Riemann invariant R \ , V t and S are 
all set to free-stream values. The other variable Rz is extrapolated from the interior flow variables. On 
outflow, V n > 0, only R\ is fixed to free stream and Rz,Vt, and S are extrapolated. Once these four 
variables are available at the boundary the four flow variables Q can be obtained. 

For these computations, the wake cut in the “C” meshes is handled implicitly. The array storage 
is shuffled so that integration in the rj -direction proceeds directly from one outer boundary to the other 
in the wake region, and from the body surface outward to the far-field boundary in the rest of the grid. 


RESULTS 


Computational results for the following five viscous transonic airfoil cases (from ref. 25) were 
compared: 

1. (Al) NACA 0012 at Moo = 0 .7 , a exp = 1 .86°, cw r = 1 .49°, Re = 9 x 10 6 , transition at 0.05 
chords. 

2. (A2) NACA 0012 at = 0 .55, a eip = 9 .86°, a corr = 8 .34°, Re = 9 x 10 6 , transition at 
0.05 chords. 

3. (A4) NACA 0012 at = 0 .7 , a » 3 .0°, Re = 9 x 10 6 , transition at 0.05 chords. 

4. (6) RAE 2822 at M = 0.725, a txp = 2.92°, = 0.729, cw = 2.31°, Re = 

6 .5 x 10 6 , transition at 0.03 chords. 

5. (10) RAE 2822 at Moo„ p - 0.750, a exp = 3.19°, M oocorr = 0.754, acoTT = 2.57°, Re = 
6.2 x 10 6 , transition at 0.03 chords. 

The angle-of-attack corrections for the NACA 0012 cases were obtained from reference 26. The 
corrections for the RAE 2822 cases were obtained from Hall (private communication from M. G. Hall, 
Royal Aircraft Establishment, Famborough, UK, 1988). In those cases, the airfoil coordinates were 
measured coordinates with a camber correction, rather than the usual design coordinates, as proposed 
by Hall. Experimental data for cases Al and A2 are from reference 26, and those for cases 6 and 10 
are from reference 27. 

Each of the two codes — ARC2D and FLOMG — were used to compute each of the five cases on 
three different grids generated by Maksymiuk, Swanson, and Hall. All grids were large C-meshes with 
leading-edge and trailing-edge clustering and very fine spacing in the normal direction near the body 
surface. The two grids (one for the NACA 0012 and one for the RAE 2822) designated “Maksymiuk” 
were generated with a hyperbolic solver (ref. 28). These grids are 385 x 65 points, with 321 points on 
the body (streamwise) and 32 points in the wake. The only clustering of points is at the leading and 
trailing edges, where the spacing is 0.001 . Normal spacing at the body is 0.00001. The outer boundary 
is at 25 chords. The grids designated “Swanson” are also 385 x 65 points, with 257 points on the body 
and 64 in the wake. The leading- and trailing-edge spacing is 0.001 on the NACA airfoil and 0.002 
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on the RAE airfoil. Normal spacing at the body is 0.000004 on the NACA airfoil and 0.00001 on the 
RAE. The distance to the outer boundary is about 20 chords. The “Hall” grids are slightly less resolved 
at 321 x 65 points, with 257 points on the body and 32 points in the wake. The outer boundary is at 
8 chords. The NACA 0012 grid has a normal spacing of 0.0000127 at the wall and clustering at the 
leading and trailing edges, where the spacing is 0.00066 and 0.00481, respectively. The RAE 2822 grid 
has a normal spacing of 0.0000555, a leading-edge spacing of 0.00064, and a trailing-edge spacing of 
0.00457. 

All of the figures showing computational results are found at the end of the paper. Figures 3-7 
show the pressure coefficients on the airfoil surfaces; results from each grid are shown separately. For 
case A1 (fig. 3), the computations show a weak shock near the leading edge. Despite the apparently 
larger amount of second-difference dissipation, the ARC2D gives a slightly sharper shock than the 
FLOMG. This behavior of the two codes is a consequence of the several differences in the numerical 
dissipation employed. In particular, there are differences in the form of the eigenvalue scaling, the 
magnitude of the constant coefficient for the fourth-difference term, and implementation of the switch- 
ing operator. In case A2 (fig. 4), there is a strong shock at about 20% chord. On the highly resolved 
Maksymiuk grid (fig. 4(a)), the ARC2D nearly captures the lambda shock structure; however, the main 
part of the shock is slightly farther aft than the experiment indicates. The FLOMG exhibits a sim- 
ple strong shock in this region. On the coarser grids (figs. 4(b) and 4(c)), both codes compute simple 
shocks, with the ARC2D locating it slightly farther aft on both grids. The same behavior is observed 
for case A4 (fig. 5). On two of the grids (figs. 5(a) and 5(c)), the FLOMG results show a smaU over- 
compression at the shock. Results for both cases 6 and 10 are similar. On the finest grid (figs. 6(a) and 
7(a)), there is again a slight difference in shock location. On the coarser grids (figs. 6(b), 6(c), 7(b), 
and 7(c)), the two codes show slight disagreement along the upper surface, as well as at the shock. The 
large error in shock location relative to experiment in case 10 has been attributed to a shortcoming in 
the turbulence model for cases with large separation zones (ref. 25). 

Figures 8-12 show skin-friction computations for each case. In the plots for case A1 (fig. 8), the 
ARC2D results show a bump at 20% chord, corresponding to the sharper shock produced there. The 
FLOMG results for cases Al, A2, and A4 on the Hall grid (figs. 8(c), 9(c), and 10(c)) have a sharp 
increase at the trailing edge. Figure 11 illustrates how important grid resolution is to the proper calcu- 
lation of skin friction. Experimental results are closely matched by both codes on the Maksymiuk and 
Swanson grids (figs. 11(a) and 11(b)); however, on the Hall grid (fig. 11(c)), both codes underpredict 
the friction coefficient. The Hall grid does not have adequate grid resolution in the normal direction 
at the body surface. Experimental results are less successfully matched in case 10 (fig. 12), but both 
codes still compute somewhat lower friction on the Hall grid. 

Selected boundary-layer profiles from cases 6 and 10 are shown in figures 13-16. Case 6 results 
are in figures 13 and 14. Both codes are in good agreement with each other and with experiment. Some 
of the dissimilarities might reflect difficulties in pinpointing the edge of the boundary layer. The results 
for case 10 (figs. 15 and 16) do not match those of the experiment, since the turbulence model precludes 
computing the post- shock separation accurately. 

Displacement thickness is plotted in figures 1 7 and 1 8. The results of FLOMG for case 6 are quite 
good, but those of the ARC2D underpredict the experiment downstream of the shock. This could be a 
result of the greater dissipation used in its computations. 
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Table 1 displays lift, drag, and convergence data, as well as CPU time in seconds needed to achieve 
a given percentage of the converged value of lift. It is evident that FLOMG, which can get a solution 
to 0.1% of the converged value of the lift coefficient in just over a minute for most cases, is much 
faster than ARC2D, which takes 5-11 min. For computing steady-state solutions, the extra work of 
an implicit method is not needed. However, the explicit scheme is significantly less efficient for un- 
steady flows than for steady ones, a result of the inapplicability of local time-stepping and the low 
effectiveness of the current form of the multigrid method. Improvements are being investigated. Some 
of ARC2D’s convergence times are unusually high because the runs were not optimized to converge 
quickly. Also, implicit wake integration operates more efficiently on the Swanson meshes, which have 
twice as many points in the wake as the Maksymiuk meshes, so there are more long vectors. Conse- 
quently, convergence times are generally shorter on the Swanson grids than on the Ames grids. It is 
interesting that FLOMG requires about the same amount of time in all cases, whereas ARC2D is quite 
variable (especially for case 10, which was difficult to converge because of slight motions of the shock 
on the finely spaced grids); moreover, ARC2D usually requires proportionately more time for the final 
stage of convergence (from 1% of Cj to 0.1% of Ci) than FLOMG. 

Some additional results obtained with FLOMG are presented in the appendix. They indicate the 
influence of artificial dissipation on several of the computations. 


CONCLUDING REMARKS 


ARC2D and FLOMG produce essentially the same results. Pressure coefficients are nearly identi- 
cal, except in some cases near shocks; some of these differences are attributable to different dissipation 
coefficients. Their agreement with experimental data is excellent, except for shock location in flows 
with large separation zones. Skin-friction coefficients, boundary-layer profiles, and displacement- 
thickness results are very similar. Satisfactory agreement with experimental data is observed for grids 
with sufficiently fine spacing near the body. 

Lift coefficients calculated by the two codes are generally within 1% of each other, and drag 
coefficients generally fall within 5% of each other. For attached cases, both codes predict lift within 
5% of the experimental value, and drag within 10%. This accuracy is not achieved for flows with 
separation, where deficiencies in the turbulence model are apparent. 

For the steady-state cases computed here, the explicit multigrid method demonstrated superior 
convergence characteristics, requiring only about 1 min of CPU time on a Cray 2, compared with 
5-10 min for the implicit method. 
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APPENDIX 

THE EFFECT OF DISSIPATION COEFFICIENTS ON FLOMG RESULTS 


In FLOMG computations, the constants for the second-difference and fourth-difference artifi- 
cial tLtotZ terms were selected because they had worked quite weU for a wtde vanety of two- 
dlltiSexlal and internal flow problems. However, a. leas, for the atrfotl cases considered 
herein, which include the difficult RAE airfoil case 10, the dissipation consumts can be ^“^y 
factor of 2 without causing serious deterioration in convergence behavtor Stnce *«erence :oc . 
in the vicinity of shock waves between some of the results obtained with FLOMG and ARC2D, several 
computations were repeated to determine the effect of the smaller dissipation factors. 

Very small changes occurred in the sharpness of the shocks for the cases considered^ As seen 
in table 2 2 maximum increase in lift was less than 3%. Except for case 10, in which the toud 
drag increased by 5 counts, the drag increased by 2 counts or less. The computer times required to 
reach the different levels of lift-coefficient convergence are essentially the same as those of the hig er 
dissipation constants This again demonstrates the robustness of the scheme. Pressure coefficients 
£ m in table 2 are panted in figures 19-21; skin-fricdon coefficients are presented 


in figures 22-24. 
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Table 1 . Lift coefficients, drag coefficients, and convergence times. 









CPU time required to 
reach designated value 

Case 

grid 

code 

Q 

Cd 

Co, 

Cd, 

of Ci convergence, sec. 
5% 1% 0.25% 0.1% 


A1 Maksymiuk 



experiment 


A2 | Maksymiuk 


ARC2D 0.253 0.0082 0.0025 


FLOMG 0.252 0.0081 0.0024 0.0056 


ARC2D 0.256 0.0088 0.0023 0.0065 


FLOMG 0.253 0.0081 


ARC2D 0.257 0.0073 


FLOMG | 0.258 0.0077 1 0.0023 1 0.0052 


0.241 0.0079 


imiitohlilEiE 


52 


25 250 


52 


380 


39 


572 678 


72 


350 415 


60 


493 556 


51 



ARC2D 0.978 0.0357 


FLOMG 0.984 0.0352 0.0314 0.0038 


ARC2D 0.994 


0.0321 0.0047 


experiment 


FLOMG 0.981 0.0353 0.0316 0.0037 


ARC2D 0.991 0.0360 


FLOMG 0.949 0.03411 0.03081 0.00331 


0.983 0.0253 


59 254 


56 260 


59 


250 


62 


335 356 


81 


310 450 


78 


350 384 


79 


ZMmiESxm 




ARC2D 0.503 0.0143 


FLOMG 


0.0085 

0.0052 

0.0089 

0.0060 

0.0086 

0.0053 

0.0088 

0.0047 


52 


23 230 


FLOMG 0.493 


6 Maksymiuk 



experiment 


10 Maksymiuk 



0.778 

0.0143 

0.772 

0.0136 

0.774 

0.0141 

0.783 

0.0138 

0.773 

0.0132 

0.782 

0.0134 

0.743 

0.0127 

0.805 

0.0292 

0.811 

0.0289 

0.800 

0.0288 

0.814 

0.0287 

0.804 

0.0279 

0.807 

0.0276 

0.743 

0.0242 


46 

112 


51 


529 661 


72 


317 385 


60 


293 318 


56 



572 695 


70 


621 745 


70 


388 502 


71 



1478 1511 


72 


379 935 


99 


126 1 474 


83 
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Table 2. Lift coefficients, drag coefficients, and convergence times 
for FLOMG results with reduced artificial dissipation. 


Case 

grid 

c t 

c rf 

C«p 

c df 

CPU time required to 
reach designated value 
of Ci convergence, sec. 

1% 

0.1% 

A2 

Maksymiuk 

0.990 

0.0352 

0.0312 

0.0039 

53 

81 

. 

Hall 

0.963 

0.0340 



62 

79 

6 

Maksymiuk 

0.779 




54 

70 


Hall 

0.794 

0.0136 

0.0082 

0.0054 

45 

74 

10 

Maksymiuk 

0.820 ! 




52 

72 


25 































x/c 


(a) Maksymiuk grid. 

Figure 3.-Pressure coefficient for NACA case Al: M 0 0 = 0.7, ot exp = 1 .86 , Q^corr — 1.49°. 
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(b) Swanson grid. 
Figure 3. -Continued. 
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(c) Hall grid. 

Figure 3.-Concluded. 
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(b) Swanson grid. 
Figure 4.-Continued. 
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(c) Hall grid. 

Figure 7.-Concluded. 
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(a) Maksymiuk grid. 



(b) Swanson grid. 



(c) Hall grid. 

Figure 13.-Boundary-layer profile at x = 0.319 for RA 
= 0.729, a eip = 2.92°, a corr = 2.31°. 





case 6: 
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(a) Maksymiuk grid. 



r 

.6 -8 


(b) Swanson grid. 



u / oo 


(c) Hall grid. 

Figure 14.-Boundary-layer profile at x = 0.95 for RAE case 6: M« 
Moo = 0.729, a exp = 2.92°, a corr = 2.31°. 





.05 


i 



(a) Maksymiuk grid. 



(b) Swanson grid. 



(c) Hall grid. 


Figure 15.-Boundary-layer profile at x — 0.75 for RAE case 10: M c 
Moo corr = 0.754, a e * p = 3.19°, a corr - 2.57°. 









(a) Maksymiuk grid. 



(b) Swanson grid. 
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(c) Hall grid. 

Figure 17.-Displacement thickness on upper surface for RAE case 6 

Moo„ p = 0.725, Moo corr = 0.729, a ea!p = 2.92°, a corr = 2.31°. 








(a) Maksymiuk grid. 

Figure 19 -Pressure coefficient for NACA case A2: M 

Q AO 


0.55, cx eX p — 9.86°, oc corr 








x/c 

(a) Maksymiuk grid. 

Figure 20.-Pressure coefficient for RAE case 6: M <*> = 0.725, Moo = .729, a ern = 

2.92°, a corr = 2.31°. 
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(a) Maksymiuk grid. 



(b) Hall grid. 

Figure 22.— Skin-friction coefficient for NACA case 2: Moo = 0.55 ,a e xp — 9.86 ,oc corr 





(a) Maksymiuk grid. 



x/c 


(b) Hall grid. 

Figure 23.-Skin-friction coefficient on upper surface for RAE case 6: M 00€XT = 0.725, M, 
0.729, a exv = 2.92°, a corr = 2.31°. 



Figure 24.-Skin-friction coefficient on upper surface for RAE Case 10 on Maksymiuk grid 

Moo... — 0.75, = 0.754, a exp — 3.19°, a corr = 2.57 . 




Report Documentation Page 

Spao* Adminitfoton 


1. Report No. 2. Government Accession No. 

NASA TM- 1028 15 

4 Title and Subtitle 

A Comparison of Two Central Difference Schemes for Solving the 
Navier-Stokes Equations 


3. Recipient's Catalog No. 

5. Report Date 

July 1990 

6. Performing Organization Code 


7 Author(s) 8 Performing Organization Report No. 

C. M. Maksymiuk, R. C. Swanson (Langley Research Center), and A-90141 
T. H. Pulliam I 1 n VA/nrlt I Init Kin 


9. Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035-1000 

1 2. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. Work Unit No. 

505-60 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 


15. Supplementary Notes 

Point of Contact: Catherine M. Maksymiuk, Ames Research Center, MS 202A-1, 
Moffett Field, CA 94035-1000 
(415)604-4737 or FTS 464-4737 


16. Abstract 

Five viscous transonic airfoil cases have been computed by two significantly different computa- 
tional fluid dynamics codes: an explicit finite-volume algorithm with multigrid, and an implicit finite- 
difference approximate-factorization method with eigenvector diagonalization. Both methods are de- 
scribed in detail, and their performance on the test cases is compared. The codes utilized the same grids, 
turbulence model, and computer to provide the truest test of the algorithms. The two approaches produce 
very similar results, which, for attached flows, also agree well with experimental results; however, the 
explicit code is considerably faster. 


17. Key Words (Suggested by Author(s)) 
Computational fluid dynamics 
Navier-Stokes equations 
Airfoils 


18. Distribution Statement 

Unclassified-Unlimited 


Subject Category-02 


19. Security Classif. (of this report) 

Unclassified 


20. Security Classif. (of this page) 

Unclassified 


21. No. of Pages 
62 


22. Price 


NASA FORM 1626 OCT 06 


For sale by the National Technical Information Service, Springfield, Virginia 22161 








